


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1972 


Stability and convergence studies ona 
numerical model of the Arctic Ocean. 


Willems, Robert Cecil. 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/16364 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 


/ (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist : Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

; | LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


STABILITY AND CONVERGENCE STUDIES 
ON A NUMERICAL MODEL OF THE ARTIC OCEAN 


Robert Cecil Willems 








NAVAL POSTGRADUATE SCHOOL 


Monterey, California 





Stability and Convergence Studies 


On a Numerical Model of the Arctic Ocean 
by 


Robert Cecil Willems 





Thesis Advisor: J.A. Galt 


March 1972 


Approved for public release; distribution unluncted. 





Stability and Convergence Studies 


On a Numerical Model of the Arctic Ocean 


by 


Robert Cecil Willems 
Lieutenant, United States Navy 
B.S., Fort Hays Kansas State College, 1965 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN OCEANOGRAPHY 


from the 


NAVAL POSTGRADUATE SCHOOL 
March 1972 








ABSTRACT 


Complete solutions to the equations of motion are possible using 
numerical schemes. The model developed by Galt (1972) was investigated 
_ for stability and convergence. A series of runs indicated that the 
scheme, the three level Adams-Bashforth method for the new value of 
vorticity and then the successive over-relaxation method for the stream 
function, was stable for transport not exceeding 100 Sverdrups and a 
time step of one-half pendulum day. Convergence was determined by 
comparisons, with known analytic solutions, of phase velocity of Rossby 
waves, and for the steady state conditions, patterns of the stream 
function and the vorticity. Improved accuracy of the model was attained 


by using a finer grid. 
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I. INTRODUCTION 


The study of the dynamics of oceanic flow has been attempted by 
uSing the equations of motion. Since solutions to those complex equa- 
tions in their complete form was not probable it’became necessary to 
Simplify them by making highly restrictive assumptions that allowed 
for some of the functional relationships to be defined as constants and 
the non-linear terms to be disregarded. Ekman (1905) developed the now 
classic boundary layer solution, known as the Ekman spiral, which 
Showed the effects of wind stress acting on a free surface. Stommel 
(1948) was able to reduce the complex equations to a balance between 
wind stress, a parameterization of friction, and Coriolis force. His 
solution dramatized the phenomenon of western intensification. 

The derivation of the vorticity equation from the vertically inte- 
grated equations of motion removed the necessity of defining the press- 
ure field at the surface and simplified the boundary conditions. Munk 
(1949) developed a solution to a simplified form of the vorticity equa- 
tion which contained no non-linear terms, was invariant with time, — 
assumed a zonal wind stress, and assumed a constant eddy coefficient. 
In both solutions the features Stommel observed in his solution were 
evidenced and, additionally, Munk (1949) observed the counter current 
region as well. 

Development of analytic solutions to the equations of motion 
necessarily demanded regular, closed boundaries, and for some forms, a 
constant depth. In each solution mentioned, an ocean of constant depth 


was assumed or, as in Ekman's solution, the ocean was considered as 





unbounded in depth. To study, analytically, ocean basins of irregular 
boundaries, with inflow and outflow, and variable depth lead to an 
impossible task. 

Early work in developing these solutions was usually attempted for 
ocean basins that were rectangular in shape. A need to investigate the 
effects of a laterally sloping boundary was recognized, as evidenced 
by the more recent work, by Munk and Carrier (1950) which was done 
on a triangular shaped basin. Studies of the effects of bottom topo- 
graphy on ocean circulation have been presented by Warren (1963) and, 
more recently, by Veronis (1966). The solutions Veronis developed 
were to the time dependent, non-advective equations with no friction or 
wind stress and fora rectangular basin. 

It was clearly indicated that with each solution a different pheno- 
menon of oceanic flow could be described using the vorticity equation 
or the equations of motion. The optimum stiuation, then, would be to 
obtain a complete solution to these equations so that instead of being 
restricted to solutions which describe only one of the phenomenon a 
complete interaction of advection, Coriolis force, friction, and wind 
stress could be studied as a function of time. 

Since a numerical solution to the equations of motion removed the 
necessity of simplification and allowed for a complete solution, various 
models have been developed. Usually the only simplifications involved 
were to specify the density field and friction as constant. More 
recently Galt (1972) described development of a numerical model of a 


constant density form of the vorticity equation. 





The non-dimensional equations that were used are represented by 


equations (1) and (2). 


Vortici ty: 


De ho eee eee g 2 i 
seh vy (224) ee tl (1) 


Stream function: 


ve(P vv) = & (2) 
Where: 
= vorticity 
h = depth 
f = Coriolis parameter 
K =  Jateral kinematic eddy coefficient 


In this 


Where: 


= SS. ss = 


The non- 


= wind stress 


= horizontal gradient operator stream function 


form the non-dimensionalized constants o and 6 were defined: 


Q 
it 


2 
Y/FL D 


oO 
iT 


K/FL¢ 


= a representative value of the stream function 


= a representative value of the Coriolis parameter 


= mean depth 


= the grid length 


dimensional wind stress t was related to t', the dimensional 


form, by equation (3). 





tT os c*( | (3) 





Where: 
ep = density 


The non-dimensional stream function was then defined: 


b= v' (FLD) = vy 


Using the non-dimensional forms of the equations the model was 
constructed. The new value of vorticity was obtained by the three level 
Adams-Bashforth differencing method. | 

The stream function was then relaxed to within acceptable limits by 
the method of successive over-relation. This scheme was applied to a 
polygonal computational molecule. 

This computational grid allowed for the specification of irregular 
boundaries and for inflow and outflow. Newtonian, lateral friction was 
parameterized in terms of the horizontal Ekman number, K/FL‘. 

For models which have been developed it has been necessary to, first, 
determine their stability and second, to determine their convergence to 
the true solution. Since complete, analytic solutions were not available, 
these models were run for conditions which were representative of 
Simplified equations to which solutions existed. Comparisons were hed 
made to assess the validity of the models. 

For some finite difference schemes it is possible to mathematically 
solve for errors generated and determine rates of convergence as well as 
Showing what bounds there must be on the spatial and temporal steps to 
ensure stability. In referring to the model Galt developed it was not 
possible to make general calculations for convergence and stability. 

As indicated by Baer and Simons (1970) a method to investigate numerical 
Stability inamodel can be effected by monitoring a conservative 


property such as energy and total vorticity. The method used to determine 
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convergence was by a comparison of the patterns of vorticity and the 
stream function with known analytical solutions. This method, for 
example, was used by Blandford (1970). 

Establishment of the characteristics of Galt's model in terms of 
stability and convergence was done in two parts; 1): For stability, a 
form of kinetic energy, total vorticity, and the maximum value of the 
stream function were monitored; 2): To check.convergence a comparison 
was made between known analytic solutions and the solutions predicted 
by the model. These comparisons used patterns of the vorticity and the 
Stream function, and of the phase velocity of Rossby waves. These 
comparisons were made with two different model configurations. 

Once the basic shape and boundary conditions were determined the 
initial input was either constant or variable depth, distribution of 
the Coriolis parameter, magnitude and direction of the wind stress, 
Stream function and vorticity patterns, and the local rate of change and 
vorticity one time step back. The parametric specification of density 
and frictional coefficients, which were assumed constant, and the 
physical dimensions of a particular configuration were determined by 
the magnitudes of the non-dimensionalizing constants, b /FLED and K/FLS. 
Additionally, for those runs when these constants were zero the 
dimensions of the model were determined by the distribution of the 
Coliolis parameter, magnitude of the wind stress, and bottom topography. 

The first configuration was a rectangle Shape that could only be 
approximated by the model due to the use of the polygonal molecules 
for the computational grid. The boundary conditions were set up so as 
to represent an enclosed basin with free slip boundaries with no 


friction in the model. The non-linear terms were not used and the wind 
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stress set to zero so that this solution could be compared to the free 
wave solutions of Veronis (1966). One run was with constant depth and 
a variable Coriolis parameter and the other run was with variable 
depth and a constant Coriolis parameter. 

The second series of runs were done using a parallelogram shape 
which had closed regular boundaries. The east-west boundaries were 
inclined sixty degrees from the horizontal, measured in a positive 
Sense. This was run with constant depth and a variable Coriolis 
parameter. The north-south boundaries were set to free slip conditions 
and the east-west boundaries were set to no slip conditions with 
friction in the interior. A zonal wind stress set in as a consine func- 
tion with maximum values at the north-south boundaries. The model did 
not accurately convert this to the curl of wind stress and, therefore, 
all successive runs were made with a direct input of the curl of the 
wind stress. These particular requirements were so that a comparison 
with a modified solution of Munk and Carrier's (1950) could be made. 

The investigation of the effect of reducing the spatial step size 
was done for both configurations by increasing the number of points 
and appropriately scaling the input parameters. Reduction of the grid 
length to one-half was accomplished so as to better define the 


velocities and western boundary layer. 
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II. STABILITY ANALYSIS 


Baer and Simons (1970) indicated that for the Adams-Bashforth 
iteration scheme, when used on a non-linear equation, it was sufficient 
to determine stability for a linear form of that equation Since the on- 
set of instability appears to be the breakdown of the linear solutions. 
Since determination of computational stability was desired it was suffi- 
cient, (Baer and Simons [1970]), to monitor the conservative properties 
which were, kinetic energy, total vorticity, and the maximum value of 
the stream function. For the stable case, Figure 1, these properties 
were bounded over the period of integration. Figure 2, the unstable 
case, indicates some of these values as growing unbounded after a certain 
integration time. The model's response for the stable situation was a 
comparatively short spin-up time, as indicated by the values of the 
properties monitored. The-kinetic energy and maximum value of the 
stream function reached a maximum value at this time and were then 
observed to oscillate throughout the remainder of the integration time. 
This response was similar to the results reported by Blandford (1970) 
for a model in which was specified a lateral viscosity and no slip 
boundary conditions. Total vorticity reached as maximum value later in 
the integration then remained constant throughout the remainder of the 
integration time. For the oe Situation the kinetic energy after 
a period of time began to grow in an exponential manner with oscilla- 
tions of the stream function increasing. The vorticity, for this 


Situation, grew steadily throughout the integration time. 


ks 





"ase) ajyqeys out “| aunbly 


AJWIL T3GOW TWNOISNSWIC-NON 


_ OF9 09S O8¢ O0r OZE Ord O9t 038 
oe 
P “a 
“ 
“ 
| ma 
pr O'l 
4 
p-OL X APLILZAOA 
_— ———" {PU LSUdWLP-UON 
UoLzoUN{ weaAzS ay? JO aNLeA es 
wNnwiXeW [| PEUOLSUSWLB-UON! |, ae ~ 
, ostoz 
ime 
. = 
f \ oO 
i T] 
‘ ae 
Tm 
/ \ ‘ = 
7 
/ ./ x OTT O'E 
fT) 
a> 
= 
NOLQOUNY weau7zs ay _ 
$O ONLeA wnwWwLxeW = 
| PEUOLSUDWLP-UON a 
= o'9t0'y 
> 
ae 
Co 
< 
oO X ABuauZ OLYOULY [LeUOLSUSWLP-UON zs 
~~ 007+0°S 
— 
~< 


14 





OzE 


082 


BUX ASL FA0Q 
| LUO LSUAWLP-UON 


§ 


"asey) ajqeysun aut °2 ounbl4 


AWIL TAGOW TWNOISNSWIG-NON 


Ore 007 O9t OZ 08 Or 
~ 
= 
> 
” 
4 
B 
¢ OL x ABusuy ILYoeULy 
7 -UuO —_ 
| EUOLSUBWLP-UON = 
a 
a 
a 
— 
— , -_ 
— 
one oa ae ae J ay) a 
| 


\ 


uoLzoun4 weau4s auy JO aNnLeA 
NWLXeW LeUOLSUaWLP-UON 





91 


£6 


ce 


OY 


2 





To more accurately determine the limits of stability various 
conditions were imposed on a constant depth model basin with basic 


configuration as shown in Figure 3. 


Wind Stress 


Figure 3. Basic Model Configuration for Stability Analysis 


For all but two runs mixed boundary conditions were used with the north 
and south boundaries being free slip. Of the other two runs, one was 
run with free slip boundary conditions, the other with no slip boundary 
conditions. 

The distribution of the wind stress over the model basin was zonal 
in nature. The direction of the wind stress was the cosine function in 
the north-south direction with the maximum stress in the westward 
direction at the south boundary and the maximum eastward stress at the 
north boundary. 

The first series of runs were made with the constant (y, /FL°D) set 
to zero which resulted in the linear form of the equations. The wind 


Stress was set at a representative value and the Coriolis parameter was 
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allowed to vary over the constant depth model basin. The physical dim- 
ensions of the basin were determined from this distribution. Table I 


Shows the stability for a range of friction, K. For the grid length 


“Feacn!, i 8 


lL = 4.6 x 10°m, fe LO ; K ranges from approximately 7.5 x 10° to 


7.8 x 10’ cmesec7!. 


The effect of varying the magnitude of the wind stress 
for this test had the effect of varying the transport, and by using 


large values for -the-wind stress the model was unstable. 


Table I. Stability Dependence on Friction and Time Step. 


K/FL® DT. KDT/FL® STABLE 
2.86 x 107° 0.5 1.43 x 10> Yes 
2.86 x 107" 0.5. 1.43 x 107° Yes 
2.86 x 10°" 1.0 2.86 x 107" Yes 
2.86 x 107" 2.0 5.72 x 104 Yes 
2.86 x 107" 3.0 8.58 x 107" ? 
2.86 x 107" 5.0 Tee 10” ? 
2.86 x 107" as Toyo. Ky No 
2.86 x 107° —-10.0 2.86 x 10° No 


Losin | 


Allowing K to remain constant at 7.5 x 10°cm sec the wind was 


varied as shown in Table II. 


Table II. Stability Dependence on Transport 





cur] i DT STABLE 
-sin(y) 0.5 No 
0.0 O25 Yes 
S107 “sinty) 0.5 . Yes 
-10° “sin(y) 0.5 Yes 
-10° lcin(y) 0.5 Yes 
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By using the non-dimensionalizing coefficient for the wind stress, 


ge 


where p is density, F a representative value of the Coriolis parameter, 


T 


X 


eal, 


X 


L the grid length, Yo in Sverdrups (10°m?/sec) was determined. The 


maximum absolute value of the non-dimensional wind stress, t, was 1.0. 


with o 


in Table II represent 


yo" grams m 


Sverdrups. 


3 


, =~F-= 107 sec” 


] 


gels, =efl OeX oliO 


mn the runs depicted 


a range of maximum transport from zero to 10° 


The effects of non-linearities on the stability of the model are 


Shown in Table III. 


Keeping K constant at 7.5 x 10 


to vary with varying wind stress. 


Inclusion of the non-linear terms by Yo /FL 
depth, y, ranged from 10° 


magnitude.of the wind stress and it ranged from approximately 2 x 10° 


Table III. 


cur | ue 


-107 sin(y) 
00 
-10° 
-10 ‘sin(y 


'sin(y) 
(y) 
-107!sin(y) 
(y) 
(y) 


] 


Sinly 


-10"' 
-107! sin(y 


sinty 


-9 2 -2 


to 2x 10 ~ msec 


/ 2 


: o/FL D was allowed 


Stability Dependence on Non-linear Terms and 


Time Step. 


2 
Yo /FL D 


2 
2 
2 
3 
3 
6 
3 


10° 
10 
10° 
oa 
10° 
10° 
10° 


2 


a=) 





DT STABLE 
1.0 Yes 
10 Yes 
10 No 
1.0 Yes 
5.0 ? 
1.0 Yes 
0.0 No 

2 


for its maximum values. 


D where D is the mean model 


to 10° Sverdrups. This also scaled the 


1S 


Similar, unstable, response 


for large transport was obtained by setting ¥o/FLAD to 10.0 which 


represented approximately 10° Sverdrups. 
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Table IV indicates the final three runs made for determination of 
computational stability. Again, using the linear form of the equations 


2 1 


and K approximately 10°cm sec. , the maximum transport for the wind 


stress used was approximately 10 Sverdrups. 


Table IV. Stability Dependence on Time Step. 


2 


Curl K/FL DT STABLE 
-107*sin(y) ale 1.0 Yes 
-107'sin(y) 5x 107" 5.0 ? 
-1074sin(y) See 10.0 No 


An attempt to compute computational stability was made by using; 


At At 
c ft 4 x IE 1.0 
( AX (ax)? 


The term, C A where C is a characteristic velocity, was not important 

: : : At 

in the series of linear runs. The term K Taxye Was computed as shown 

in Table I with the results of a somewhat larger bounds for At. For 

the transport range of 10 to 100 Sverdrups and a coefficient of friction 


sae 4 ] 


on the order of 10°cm°sec 


a time step of 10° seconds ensured stability. 
This represented a one-half pendulum day. The effects of the different 
boundary conditions mentioned previously was to change only the model 


Spin up time. This observation was similar to that reported by Blandford 


(1970). 
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III. ROSSBY WAVE PROPAGATION 


Convergence of the model to a steady state analytic solution was 
determined by allowing the model to generate time dependent, free wave 
solutions and allowing the transients to die out leaving only the steady 
Solution. These Rossby waves can occur in an ocean basin and are of 
two types. The classical Rossby wave can exist in a constant depth 
ocean due to the change in the Coriolis parameter with latitude. Con- 
Sidering the situation where there is no friction and no wind stress, 
with the Coriolis parameter increasing from south to north, a Rossby 
wave will tend to propagate along lines of a constant Coriolis parameter. 
For this case, then, the wave would propagate in an east to west direction. 

The second type of Rossby wave that can exist in an ocean basin is 
the topographic Rossby wave. Consider, now, a constant Coriolis para- 
meter but with a simply varying depth which increases from east to 
west. The topographic Rossby wave will propagate as before except now 
along lines of constant depth in an east-west direction. For a simply 
varying depth, increasing from east to west the topographic Rossby 
waves will propagate from south to north across the ocean basin. These 


waves can be described by solutions to the potential vorticity equation, 


e 


ae 4 vey (f)-; —E = vorticity (4) 


‘) 


wv 


ctr 


Considering that potential vorticity is conserved following a parti- 
cle of water’, equation (4) states that for free wave solutions the second 
term is zero, tat is, f/h is constant. Using equation (4), Veronis 


(1966) developed a time dependent, harmonic, constant depth solution 
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for the stream function, w (x,y,t). 


y= 


0 


sin( Ax + wt) sinox sinyy (Veronis[1966]) (5) 


Equation (5) describes sinusoidal waves propagating across a basin in 


the -x direction enclosed in a sine envelope. For the first mode in the 


x and y directions 
A 


a 


M 


the following constants were defined: 


= p/2 
= n/L 
= 7/M 
= p/2far + y77* 
= constant = ar 


~ basin dimension in x direction 


= basin dimension in y direction 


For use in initializing the model non-dimensionalization of equation (5) 


was necessary. By 


non-dimensionalizing L, M and A with a grid length, 


AL, and w with F, the following non-dimensional constants were defined: 


The phase velocity 


; 2 tah 
= [a “ ty “J 


= (n/L) AL 
= (71/0.866L) AL 
= BAL/2a F 


grid length 
= basin dimension 


= a representative value for the Coriolis parameter. 


for this solution is w/A which means that the non- 


dimensional phase velocity is w'/a' = BAL /2x° F, or gAL/2Flo - Y al 


- A second solution developed by Veronis (1966) was for a linearized 


form of equation (4), but with a constant Coriolis parameter and with 


2] 





variable depth that was assumed to be h = hee. This solution, 


equation (6), described sinusoidal wave patterns: 


v= Y ekx/2 sin(ay + wt) sinax sinyy (6) 


propagating in the positive y direction which are enclosed in a sine 
envelope. The envelope is increasing in an exponential manner in the 


x direction. Again, for the first mode, the dimensional parameters were 


defined. 
A = Fk/2w 
0 = afl 
yY = 1/M 1 
w = Fk/(ALo? + y2] + k2)2 
f = constant Coriolis parameter 
L = basin dimension in x direction 
M = basin dimension in y direction 


Using AL to non-dimensionalize L, M, and A, andw with F the following 


parameters were defined: 


kK = KAL 1 
» = kY(4[a 2 +y 2] +k 22 
“re (n/L) AL 

¥ = Gr/0-866L) AL 


The non-dimensional phase velocity,w /A , was obtained: 
I | | J | I i 
1 NDE Ui ce = eae ene 
The rectangle shaped basin which was desired could not be exactly 
Specified as depicted in Figure 4. This model basin shape had sawtoothed 


east-west boundaries due to the computational grid being made up of 


Ze 








Ate ri 


| J 


1 | ‘ 
|| a et 








polygonal molecules. Scaling the model basin to physical dimensions was 
done in the following manner. For the constant depth configuration the 


model was scaled on the distribution of the Coriolis parameter. 


ae) 
| oF _ AF 
AM FS eee ie -_ 
I a) oe a ae cy ey 
AF = 0.75 x 10 ‘sec 
a = 60° 
using a representative value for B = 
as 1.8 x 107! > aece. 
AM = = 4.63 x 10° meters 
(9)(1.8x10 —°) 
AM 


AN = = 5.35 x 10° meters 


sin 60° 


The configuration for which the Coriolis parameter was constant and 
the depth varied as in equation (6) determined the physical dimensions 


of the model in the following manner: 





- -kx 
h = hoe 
h 
In iB 
xX = =naAN , n= 9,19 
k k= 2.2x10'n | 
ho = 3.0 x 10°m 
h = 1.0 x 10°m 


mie 5856. On 


AM = sin 60°(AN) = 4.83 x 10°m 
Since free wave solutions were desired, all of the boundary condi- 

tions were set at free slip. The constants K/FLS and v /FLOD were set 

to zero so that the model had no friction and so the linear equations 


could be solved. Since the wind stress was also zero for these runs 


initialization of the stream function and corticity patterns was 
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necessary. To accomplish this, the non-dimensionalized solutions 
developed earlier from equations (5) and (6) were solved at the non- 
dimensional time 5 and the resulting patterns used as input to initialize 
the model. 

In general, the model propagated the two types of Rossby waves across 
the basin with the kinetic energy and total vorticity remaining essen- 
tially constant. The fluctuations in kinetic energy that were observed 
had their maximum value when the wave, depicted in Figure 4, reached 
the center of the basin. This was expected since the waves were pro- 
Ppagating in a sine envelope which reached its maximum at the center of 
the basin. The influence of the boundaries on the waves appeared to 
cause them to disperse laterally. This observed influence quickly 


disappeared as the waves moved to the interior of the model 


<——__——._ Direction of Movement 


Figure 4. Example of Rossby Wave Propagation. 
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Reducing the grid length to one-half its former value by doubling the 
number of points on a side did not change these observed characteristics. 
Subsequent analysis indicated that these waves were more accurately 
defined with the decreased grid length. 

The phase velocity for both series of runs was determined by using 
Figures 5 through 8. The first comparisons with analytic solutions are 
shown in Table V. This was for the configuration denoting a basin of 


constant depth and with variable Coriolis parameter. 


Table V. Comparison of Rossby Wave Velocity. 


‘NUMBER OF GRID GRID LENGTHS/10 NON-DIMENSIONAL PHASE Vel 
POINTS/SIDE TIME STEPS Model Analytic Solution 
10 Tez .12 169 
20 2.4 24 356 


The percent error for ten grid points was 29 % and that for twenty grid 
points was 17 %. 

The configuration of a basin with an exponentially varying depth 
and with a constant Coriolis parameter was analyzed in a similar manner 
for topographic Rossby waves. JTable VI indicates the comparison made 


with the analytic solution. 


Table VI. Comparison of Topographic Rossby Wave Velocity. 


NUMBER OF GRID GRID LENGTHS/10 NON-DIMENSIONAL PHASE Vel 
POINTS/SIDE TIME STEPS Model Analytic Solution 
10 1.2 ma 2 212 


20 320 35 »445 
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The same trend was observed in that the percent error for ten grid 
points was 33 % and that for twenty grid points was 21 %. 

This series of runs established a very important trend for the 
model. Increasing the number of computational points increases the 
accuracy of the model for the free wave solutions. In approximate 
terms, doubling the number of computational points very nearly halves 
the percent: error as compared to the known analytic solutions. This 
determination was subject to graphical interpretation which introduced 
possible errors. In an attempt to clarify the picture, the figures 
were normalized with respect to the sinax envelope but no significant 


improvement resulted. 


30 





IV. THE STEADY STATE SOLUTION 


The assumption that steady state conditions exist in an ocean basin 
was made by Munk and Carrier (1950) when they derived a beta plane solu- 
tion on a triangular basin of constant depth. This solution to the 
vorticity equation balanced frictional effects, Coriolis force, and an 
assumed wind stress distribution. To make this compatable with the 
configuration of the numerical model, Figure 9, the dependence of the 
basin width in the x direction on y was removed so as to make it a 
constant width. The solution was then corrected for the boundary 


conditions which were no flow normal to the boundaries. This solution, 


Wind Stress 
Figure 9. Basic Model Configuration 


in non-dimensional form is shown in equation (4); 


y = siny[-1.158{(I-e)cos(*S pa + .524) -esin(42 pa)}e(-Pa/?) 


ae ae) (7) 


3] 





with the following parameters defined: 


e = y/P 

Cm (ess tante) </> 

6 = 60° 

@ = 1-1/3(o- ,) tang 

a5 latitude at which y = 0 
¢ =  Jatitude at which y =n 
a = x/y - (y/y) tane 

yo = width of the basin 


The modified solution is shown in Figures 10 and 1]. These patterns 
were then compared to the steady state reached by the model. To 
parallel the model with this solution the boundary conditions were 
Set at free slip conditions for the north and south boundaries and no 
Slip for the east and west boundaries. The constant, K/FL°. Cleon 
the wind stress, and the distribution of the Coriolis parameter were 
Scaled with respect to each other so as to represent a constant depth 


Sn in the x 


basin with approximate physical dimensions of 5 x 10 
direction and 4.5 x 107m in the y direction, with the curl of heard 
Stress being distributed as 10° *sin(y) over the basin. The representa- 
tive value for friction, K, using these values, was 10°cm@/sec. 

As indicated by Figure 12 this configuration reached steady state 
at about 1 x 107 seconds. This time is when the maximum value of the 
Stream function was essentially at its steady state value with the vor- 
ticity and kinetic energy within 95 % of their steady state values. To 
determine the convergence of this to the analytic solution, patterns of the 


vorticity and of the stream function were drawn, Figures 13 and 14. The 


difference between these steady state patterns that the model predicted 
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and that of Munk and Carrier could probably be attributed to the differ- 
ence between the physical conditions which brought the two solutions to 
a steady state. A numerical model which Gates (1970) used to solve a 
Slightly different form of the equations of motion produced similar steady 
state patterns. He attributed this to Rossby wave reflections from the 
boundaries with the most important reflection being the one from the 
waves which propagated’ from east to west. These waves diminished in 
magnitude with each crossing as well as diminishing in magnitude as they 
propagated through the model. The basic differences, not considering 
the finite difference methods, between Gates' model and Galt's model 
were: all boundaries no slip as opposed to mixed free slip/no slip 
boundaries, and free surface as opposed to a rigid lid. These differ- 
ences did not cause the two models to vary to any great extent other 
than a longer spin up time and a phase shift in Galt's model. The most 
important parallel is that of the type of friction which was lateral 
and did not include any bE ile eons 

Blandford (1970), after considering different combinations of 
boundary conditions, was able to show by using either free slip or no 
Slip lateral boundaries that bottom friction, expressed as an exponen- 
tial boundary layer, dissipated more energy than did lateral friction. 
By including this boundary layer in the solution Blandford noted that 
the Rossby waves which were generated with model spin up tended to be 
damped out so that they did not dominate the interior flow. The boundary 
layer also caused a reduction in the time at which the model reached 
Steady state conditions. 

The response of this model, early in the integration, was to develop 
maximum transport near the center of the basin. This cell soon propagated 


to the western edge where it dominated the flow as the boundary layer. 
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After the boundary layer was established, Rossby waves began to initial- 
ize in the southeast corner and propagate to near the western boundary. 
These inertial waves remained near the southern boundary and became 
established in the flow as the weaker counter current region. Gates 
reported observations on these waves obtaining the phase velocity of 

~ 5.30 meters/sec for their westward component. These waves, depicted in 
Figures 15 and 16, had as their western component of the phase velocity 
5.34 meters/sec for the larger grid and 4.83 mer ereeee for the smaller 
grid. 

To better define the phenomenon observed in Galt's model, a reduc- 
tion of the spatial step size was attempted. Already determined, in 
Section II, was the improved accuracy of the solution with a smaller 
grid length. The response of the model for the stream function for 
different grid lengths, AL, was observed in Figure 17 as being a higher 
maximum value for (1/2)AL. This difference was not observed in the 
comparison of kinetic energy, Figure 18. The difference, in the stream 
function, between 1AL and (1/2)4L was attributed to the influence of 
transients for the larger grid length. A profile of the stream func- 
tion for the two grid lengths, Figure 19, indicated that for the larger 
grid length, 1A4L, the solution oscillates eastward of the counter 
current region. The smaller grid length, (1/2) 4L, had a well defined 
boundary current and the weaker counter current region. The transport 
decreased rapidly eastward of this region with relatively little 
oscillation. Comparing this profile with the profile of the modified 
Munk and Carrier solution, Figure 20, it was observed that the model 
did in fact compare favorably with the analytic solution in its pre- 


diction of the western boundary current and the decrease in transport 
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eastward of the counter current region. The contoured values of the 
stream function, Figure 21, supported this observation. 

The steady state solution generated by the model could be used to 
acequately indicate oceanic transport. The accuracy of this solution 
was dependent on grid length, the more accurate solution being obtained 


~ using the smaller grid length. 


46 





™~ 
“== iam wat | wae 


-.005 





Figure 21. Dtstribution of the Non-dimensional Stream 
Function for 1/2AaL. 
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Y. CONCLUSIONS 


The model was sufficiently stable for the lengths of integration 
time used. The model had to be scaled to reflect a transport from zero 
to 100 Sverdrups with a time step which represented a one-half pendulum 
day to ensure stability. To reach steady state conditions it was neces- 
sary to depend on this factor so that the transients would be damped out. 

The curl operation was not accurate enough to produce the correct 
results when operating on the wind stress. 

Rossby waves, either topographic or inertial, were satisfactorily 
described by the model. The phase velocity was more accurately defined 
with a refined grid. Although these results were dependent on graphical 
interpretation this observed trend persisted through the following 
series of runs in Section III. The errors were significant enough to 
warrant future investigations. 

The steady state solution generated by the model compared favorably 
with the modified Munk and Carrier solution. To reach this solution a 
long integration time was necessary. The spin up time was influenced 
by mixed boundary conditions which caused a longer time to reach steady 
State conditions. The transient solution the model developed was 
Similar to that which was reported by Gates. This dominating influence 
could be more accurately described by including some form of bottom 
friction. These transients became imbedded in the flow and were signi- 
ficant in establishing the counter current region at steady state 
conditions. The oscillations observed in the stream function profile 
for larger grid length were due to these transients dominating the 


Steady state flow. 
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APPENDIX THE FINITE DIFFERENCE EQUATIONS 


The specific details of the derivation of the finite difference 
equations that are given below can be found in Galt (1972). This deriva- 
tion incompasses a combination of work done by Winslow (1967), Sadourny, 
Arakawa, and Mintz (1967), and Williamson (1969). By considering a 
computational molecule, Figure. 22, the solution for each term in equation 
(8) was transformed from a surface integration to a line integration by 


using Stokes theorem. 


cee aeons tate 2 (=) 
are h veV,, ( h SY ah thy |e (8) 
(1) (2) (3) 
The value of term (1) at point L was considered as an average of the 
potential vorticity over the molecule times the value of the stream 


FuNCctION: 


I I 
potential vorticity (py!) = on - F 
h 


ioey (: r+ r) 2 gt HL 4) (ey ' | . 
ee yer) ee NY um _—— = 3 
373 
ea eee oe 
(vio . -N _ yl- “et (Py + PY \4 
a re 
L-1 L-N 
(yi-t -| yon) (20 + Py ), 
laa 
L+N-] L-] 
(4-1 4b Lt (Py + PY ) 


.?. jLtN- -) (vt + ae 
" Sou 


: (eae jis) fovi*? + pyltn ) 
= a 
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The value of term (2) at point L was considered to be the average 


value over the inner shaded area of Figure 22. 


éucE = sft fh pL-NtL + gL-N Ft pho f ml 


pLtN-1 LAN _ cet 


Term (3) considered the average values of magnitude and direction of +t 


along each side of the molecule, then stored as a constant array since 


it waS invariant with time. 
oa = + - 
eit. 2 hi ttl . a N+] x al ] in al N+] : 
Hh sal? “LAT ,L-N+] 2 
L-N+] -N L-N+] L-N 
L+] ] T Q 7 0 
“|e ie vt) 
a el an : oy 
gol - - ete bulls ] 
> “Leip cos 
h nit 
L+N-] gl tN- ] aL tN 
2 (: L+N- 3) fn “), - ym) aL tN- 7 
h 
2 Ltn ,Ltl 2 


The successive over-relaxation procedure was used for equation (2) 




















where the residual was computed by equation (9) 
Res = of 9 v) -€ (9) 


Res was determined by an average value of V( r VY w) over the inner 


hexagon, Figure 22. 


50 





ctor)? (sate | - c 
pLtn-1 rt ne Ltn A nt 


The amount of over-relaxation was determined by equation (10). 


yt = yl + R(dy“) R = 1.48 (1.0) 


ie = Res 2 
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